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The technique now regularly in use for the generation of orbital elements 
for the TELSTAR communications satellites is described using angle-only 
and/ or angle-range data versus time as input information. It is found that 
secular perturbation considerations are sufficient to permit trajectory pre- 
diction with pointing errors of about 0.05° over 100 orbits from a single set 
of elements. Modified orbital elements are chosen as the orbit description, 
since they explicitly express the secular rates and thereby simplify and re- 
duce the cost of drive tape generations for the Andover ground station. The 
rates are derived both from perturbation theory and from direct measurement. 
The size and, shape accuracy of the predicted orbit ellipse is improved by 
statistical means using trajectory data from a number of passes over a par- 
ticular ground station. This permits better round-the-world predictions. 
Finally the computer -operator ensemble presently used for generating the 
orbits of the TELSTA R satellites is described. 
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I. INTRODUCTION 

The calculation of elements for orbits similar to those of recent near- 
earth satellites such as the first and second TELSTAR communication 
satellites, Echo I and II, the Tiros series and the Syncom satellites is ac- 
complished in this analysis by considering the secular effects of the 
earth's gravitational potential upon the satellite's orbit. Here, the earth 
is taken as an oblate spheroid whose field is independent of longitude 
and is symmetrical about the equatorial plane. In the analysis, the gravi- 
tational effects of all bodies other than the earth are ignored, as are 
atmospheric drag effects and solar light pressure. 1 For orbits at distances 
in the thousand-mile range, luni-solar attraction can change the satellite 
height by no more than a few hundred feet. This can alter ground station 
pointing angles by less than 0.01°, which is below the resolution and 
tracking accuracy of the horn antenna* at Andover, Maine which sup- 
plies the bulk of the basic pointing data for Bell Telephone Laboratories 
use. The principal planetary perturbation is from Venus, which peaks at 
3 to 4 orders of magnitude less than the luni-solar effects. Air drag ef- 
fects at 1,000 miles are at least an order of magnitude below the solar 
gravity perturbations. 

In the technique to be described, modified orbital elements (osculating 
elements plus explicit secular perturbations) are derived essentially from 
three sets of observations of the satellite. If only angle information is 
available, a modified form of Gauss' method is used to ascertain the 
range, f If some range information is available, this is added to the analy- 

* This antenna 1ms a tracking jitter of about 0.005°. Its beam can resolve the 
position of a satellite to within about 0.02°. 

t It is recognized that the range to an orbiting celestial body cannot be ascer- 
tained if the three measured sight lines from observer to the body are coplanar. 
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sis. If range data for all three observations are on hand, the modified 
Gauss method is bypassed and elements are calculated from geometry 
and perturbation theory. 

The degree of goodness of elements so obtained is ascertained by com- 
paring resulting predictions of the satellite trajectory with actual data 
from the same pass as well as with data from passes prior to or after 
the one being used to generate the elements. It has been this author's 
experience that the fit of the predicted trajectory to the orbit from which 
it came is generally within to 0.02°, but that the anomalistic period 
determination from single-pass data is not sufficiently accurate to hold 
the errors under 0.1° over more than a few orbits. However, this difficulty 
is solved by calculating two or more sets of elements from passes separ- 
ated by several days and combining results to ascertain the period and 
other secular perturbations more exactly. Table I shows the resulting 
improvements under such a mode of operation for the second TELSTAR 
satellite. 

It is to be noted that such neglected perturbations as those from the 
higher-order terms in the earth's potential function or solar radiation 
pressure may be operationally introduced by not only calculating the 
period from two passes separated in time, but also computing changes 
in the other orbital element values. This has been done in the case of 
Echo IT to enable predictions to within 0.1°* over a time of about ten 
days. 

The first portion of this paper presents techniques for developing 
orbital elements suitable for predicting future satellite pointing angles 
for the ground station collecting the original data. Later portions of 
this paper will deal with techniques aimed at producing a good round- 
the-world orbit fit and element rate improvement. Finally, a computer- 
operator ensemble used to generate modified orbital elements for the 
TELSTAR satellites will be described. 

II. EFFECT OF DATA ON ANALYSIS 

While the use of modern computers makes it entirely possible to 
record and process hundreds of pointing angles for each pass of a satellite, 

In this case such a situation cannot long prevail because of the rapid motions of 
near-earth satellites and their orbits, as well as the rotation of the ground station 
on the spinning earth. Even when the orbital inclination equals the geocentric 
latitude of the station, only for portions of certain passes is the station in the 
orbit plane. The three-point orbit analysis presented here is therefore not greatly 
restrictive. 

* This is the resolution level of the equipment which recorded the Echo II 
pointing angles. 
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Table I — Improvement Obtained by Using Several 
Sets of Elements 





Initial Pass 


Pointing Errors (Angular Offset) 


Period Determined by 


10 Orbits 
Later 


100 Orbits 
Later 


1000 Orbits 
Later 


Data from 1 pass 

Data from 2 passes separated by 

1 day 
Data separated by 1 week 


0.02 
0.02 

0.02 


0.5° 
0.03° 

0.02 


5° 

0.3° 

0.05° 


over 10° 
1° 

0.3° 



it is the purpose of this paper to describe techniques which can operate 
with a sparsity of data. The reduction of the quantity of data to be 
handled, even by automatic means, is of increasing importance both 
from a time and cost basis as more and more communications satellites 
and satellite systems become available. Minimizing the quantity of data 
is of paramount importance in implementation of small computers at 
tracking sites which often have very limited storage capacities, long ac- 
cess times, and slow machine cycles. Thus we shall be concerned with 
generating modified orbital elements from passes having as few as 6 or 
12 data points available. Under these conditions, no statistical processing 
of data prior to their use in orbit generation is advisable or in many 
cases even sensibly possible. Rather, the procedure is to first generate 
elements from sets of three pointing angles (and ranges if available) on 
a pass, and then compare trajectories produced from these elements to 
the entire original set of data. By appropriate selection of triads, data 
points unsuitable for orbit prediction may be easily discovered and 
discarded, even though it may be difficult to preselect data on a purely 
statistical basis. For use in long-range predictions, orbital element rates 
(orbital period and other secular perturbations) are determined by com- 
parisons between elements of different epochs, and these rates are used 
in place of the theoretical perturbations calculated in generations from 
a single satellite pass. 

The tracking data used to determine the orbits of the TELSTAR 
satellites come from a single ground station. Any increase in station noise 
or malfunction of trajectory equipment can therefore directly affect the 
generated elements. Because of this, it was decided not to depend upon 
an uninterrupted or near uninterrupted flow of "valid" data from the 
ground station in any analysis but to maintain flexibility by developing 
elements which would predict trajectories with reasonable accuracy over 
a number of weeks. Concentration upon the establishment of a con- 
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sistent set of secular rates derived from elements generated at various 
epochs seemed a logical approach to this end. These rates are virtually 
constant for orbiting particles like TELSTAR which lose energy at an 
extremely low rate and whose oscillatory departures from the secular 
rates are not serious. 

III. THE MODIFIED ORBITAL ELEMENTS 

Trajectory generation is simplified by choosing a set of orbital elements 
in which the only perturbations are secular ones and these arc explicitly 
expressed. This in essence reduces the trajectory prediction process to 
one of simple geometry. Specifically, the perturbations considered are 
those which produce nodal regression, perigee advance, and period modi- 
fications. They have no secular effect on the orbit inclination, eccentric- 
ity, or radius of perigee. Cyclic perturbations are ignored. A further 
simplication results by quoting the elements for an epoch that is the 
time when the satellite arrives at perigee. 

No satellite moves precisely in an ellipse about the earth, and yet the 
elements to follow will imply plane elliptical motion. It may be helpful 
to consider the satellite to be moving along an ellipse in a plane which 
itself regresses westward due to the earth's oblateness. Superimposed 
upon this motion is a rotation of the major axis of the orbit ellipse about 
the earth's center located at one focus of the ellipse. Under these condi- 
tions, the orbital elements merely state the instantaneous position of 
the ellipse in space at the epoch quoted. The modified orbital elements 2 
(commonly referred to as MOES) to be used in this analysis are: 

Epoch — the year, month, day, hour and decimal minute at which the 
elements are quoted. This epoch corresponds to the time that the satel- 
lite arrives at perigee with only secular perturbations considered. 

Inclination — tilt of the orbit plane in degrees referenced to the earth's 
equator and measured counterclockwise from the equator to the orbit 
plane as seen from a point above the ascending node. 

West longitude of the ascending node — the west earth longitude of the 
ascending node of the orbit at epoch in degrees. 

Prime sweep interval (PSI) — the apparent orbit regression expressed 
as a time. Specifically, it is the time required for a point on the earth's 
equator to pass under the ascending node on two successive occasions. 
It includes, therefore, the true regression referenced to an inertial frame, 
plus the apparent regression of almost one degree per day due to the 
earth's revolution about the sun. Often the prime sweep interval is 
expressed as one mean solar day (of 1440 minutes) minus M , where M 
is simply the number of minutes added to one day to yield the PSI. 



608 THE BELL SYSTEM TECHNICAL JOURNAL, APRIL 1965 

Argument of perigee — the angle in degrees measured at the earth's 
center in the direction of the satellite's travel and lying between a radius 
vector to the ascending node position and one to perigee. 

Perigee motion — the secular motion of perigee in degrees along the 
orbit per anomalistic period. It is positive if perigee moves in the direc- 
tion of the satellite's motion. 

Anomalistic period — the time in minutes required for the secularly- 
perturbed satellite to move from perigee to the next occurring perigee. 
Generally, during this time the perigee is also moving. 

Change in period — the change in the anomalistic period in minutes 
per anomalistic period. 

Eccentricity — the eccentricity of the orbit ellipse. A constant under 
the assumed conditions of only secular perturbation and no orbit decay. 

Radius of perigee — geocentric distance of the satellite at perigee in 
statute miles. 

IV. TRAJECTORY PREDICTION 

Since each TELSTAR satellite was launched, it has been a system 
requirement to predict the position as seen from Andover, Maine on a 
rather continuous basis. While the demand for predictions of the first 
TELSTAR satellite trajectories has decreased because that satellite is 
no longer active, the second TELSTAR satellite requirements have con- 
tinued. Predictions currently are made for every pass of the second TEL- 
STAR communications satellite according to the methods described 
herein. As a result of this type of operation it became imperative to de- 
velop a computational technique which could rapidly and economically 
generate trajectories on magnetic tapes capable of pointing an antenna 
with an accuracy of about 0.05°. The use of modified orbital elements 
for this purpose permitted such predictions by purely geometric means 
and without recourse to any integration procedures, numerical or other- 
wise. The method to be outlined here has been embodied into a program 
in regular use on an IBM 7094 computer to produce drive tapes for the 
TELSTAR satellites. The program execution times for a typical hour- 
long pass in which pointing angles are given every minute is under 0.002 
hour and costs about $1 to compute. The total over-all machine cost of 
an hour-long drive tape is about $2, which includes the use of some 
auxiliary equipment. 

The method of predicting pointing angles from a ground station will 
now be described, using the modified elements as the orbit description. 
The general procedure will be outlined without explicitly stating all of 
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the geometry, which involves mathematics no more sophisticated than 
spherical trigonometry. 3 

The elements specify the west longitude position of the ascending node 
and the argument of perigee (P) along with the secular rates so that, 
for each instant of time, the position of perigee is established as shown 
in Fig. 1. 

Since the anomalistic period is given, the true anomaly V of the satel- 
lite must always be known. This, along with the orbit inclination i, de- 
scribes the subsatellite latitude and longitude as functions of time. Since 
the radius of perigee, the eccentricity, and the satellite's true anomaly 
are available, the distance of the satellite above the earth's center is 
easily determined. 

If the latitude, longitude, and height above mean sea level of the 
ground station are known, it becomes a matter of geometry to compute 
the azimuth, elevation, and slant range of the satellite from that station. 
Fig. 2 illustrates this. The great circle arc between SSP and the ground 
station can be computed to yield the central angle g; therefore the slant 
range may be found by simple triangulation, since r and R are known 
(see Appendix B for R calculations). 

A plane tangent to a spherical earth is now inserted at the ground 
station. Azimuth in this plane is determined by solving the spherical 




Fig. 1 — Geometry for determining the subsatellite point. 
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Fig. 2 — Satellite azimuth-elevation details. 

triangle N, SSP, ST A (station). Elevation is computed from the tangent 
plane perpendicular to R. Predicted azimuth and elevation referred to 
the oblate earth are derived from the azimuth and elevation related to 
a spherical earth by tipping the tangent plane north by an appropriate 
amount depending upon the latitude of the ground station. * 

V. DETERMINING SATELLITE RANGE FROM ANGLE DATA 

5.1 A Brief History 4 

The mathematical determination of orbits for celestial bodies began 
with Newton in 1680, when he calculated a parabolic orbit for a comet 
that appeared in that year. This was perhaps made feasible by Doerfel, 
who suggested that the sun lies at the focus of a cometary orbit. Euler 
in 1744 added to Newton's work the technique of mathematically relat- 
ing time to position along a parabolic, heliocentric trajectory without 
prior knowledge of the orbital elements. Around 1770, Lambert added 
many geometrical approaches to the determination of parabolic orbits 
in an attempt to reduce the problem to one unknown. Shortly there- 

* The amount of northward tilt of this tangent plane is equal to the difference 
between the geodetic and geocentric latitude of the station and is approximately 
11 '35.6635" sin 2<p - 1.1731" sin 4<p + 0.0026" sin 6^, where <p = the geodetic latitude 
of the station. See Ref . 5, p. 480. 
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after, in 1783, LaGrange added the important assumption that all ob- 
served positions of an orbiting celestial body (comet, planet, or natural 
satellite) lie approximately in a single plane. He also introduced a per- 
turbation analysis to determine departures from the assumed plane. 
LaPlace in 1780 showed a rigorous solution 6 for an orbit formed from 
any conic section about the sun, using three positions of the celestial 
body involved as well as their first and second derivatives.* This method, 
though elegant, has not come into general use because of the difficulty 
in determining the position derivatives to sufficient accuracies. 

Olbers in 1797, analytically building upon Lambert's geometry, de- 
veloped a straightforward method for the calculation of parabolic orbits. 
The method, though not as rigorous as that of LaPlace, is sufficiently 
accurate and convenient to be used with little change to the present day. 

To Gauss belongs the credit of developing the first coordinated ap- 
proach and practical solution for the case of elliptical orbits. 7 On Janu- 
ary 1, 1801, in searching for possible planets existing in the gap between 
Mars and Jupiter, the astronomer Piazzi discovered a minor planet later 
named Ceres. It was observed for a month before it was lost in the sun's 
glare. A crude projection of its motion, as well as that of the sun, in- 
dicated it would not be visible again until about October of that year. 
Since Ceres shines at only eighth magnitude and appears starlike through 
any telescope, the search for this minor planet would have been extremely 
difficult in those days without a good prediction of its position. Within 
months, Gauss developed a method for determining an elliptical orbit 
from a sparse bit of "angle-only" data. In his method he reduced the 
problem to that of two unknowns instead of one for convenience. A re- 
sulting eighth-order equation involving range and similar to that pro- 
duced earlier by LaGrange was solved by Gauss by trigonometric means. 
The Gaussian method is essentially one of successive approximations 
which rapidly converge to the true range of a celestial body by ap- 
propriate iterations. With the exception of the trigonometric form of 
the solution that has become unnecessary in the days of high-speed digital 
computers, the method of Gauss has survived with little modification to 
the present day. The only significant computational change in the tech- 
nique occurred in 1851 when Encke recast the method in Cartesian form 
and concerned himself with rectilinear coordinates rather than the six 
orbital elements of Gauss. Even now, when LaGrange's planetary equa- 
tions may be numerically integrated 91011 and data selected to obtain 
orbital elements in many chosen forms, the simple form of Gauss' 
solution is amenable to the addition of closed-form perturbations in an 

* For a modern summary of LaPlace's method see Ref . 8, p. 168, or Ref. 9, p. 40. 
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iterative and rapidly converging mode of machine computation which 
is most conservative of time. 

5.2 General Procedure 

The adaptation of Gauss' method in the present paper will make use 
of the following data and assumptions: 

(t) At least three azimuth-elevation predictions of the satellite cor- 
related to time are available. 

(it) The ground station position is known (in terms of latitude, longi- 
tude, and height above mean sea level) . 

(Hi) The ground station is not in the orbit plane for any of the 
data used. 

(iv) The instantaneous orbit plane contains the center of mass of the 
earth. 

The computational procedures to follow are now summarized. The 
first approximation of the range to the satellite is calculated using the 
above information, a refinement made to this first range estimate, and an 
initial set of orbital elements developed. Next, perturbations are calcu- 
lated from the elements, the range is again refined, and so forth until 
no subsequent change in computed range values is observed to a chosen 
number of places or until this loop has been traversed a specific number 
of times. A final calculation of elements is then made and the resulting 
computed satellite trajectory compared to the measured one as a final 
test of the validity of the elements. Fig. 3 illustrates this procedure by 
a flow chart. 

5.3 Synthesizing the Range 

We begin the detailed synthesis of the range* by referring to Fig. 4, 
which illustrates an orbit plane about the earth along with three positions 
of the satellite assumed in that plane. A geocentric right-handed Carte- 
sian coordinate system is also shown, with the z axis coinciding with the 
earth's spin axis, the x axis pointing toward the vernal equinox f and the 
y axis orthogonal to these. In this system, which does not rotate with 
the earth, the orbit plane may be expressed as ' 

ax + by + cz = 0. v . (1) 

Since the three observed positions of the satellite are initially assumed to 
lie in this plane, each of these points must satisfy (1) and we may write 

* See also Ch. 5 of Ref. 9, Ref. 12, p. 97, and Ref. 8, p. 175. For purely circular 
orbits see Ref. 13, p. 141. 
t See p. 37 of Ref. 14. 
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Fig. 3 — Flow chart for MOE generation. 




Fig. 4 — The orbit and Gaussian triangles. 
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«.(-! + byx + czi = 
ax 2 + by 2 + cz 2 = 
ax 3 + by z + cz 3 = 
from which the determinant of the coefficients becomes 



A = 



Xi 2/1 


2l 


^2 2/2 


22 


X 3 2/3 


2a 



= 



(2) 
(3) 
(4) 



(5) 



permitting (2), (3), and (4) to be rewritten as 

(2/223 - ztVa)xi - (yiza - ZiVz)xi + (yiZi - znj 2 )x 3 = (6) 

{z 2 x 3 — x 2 z 3 )yi - faft — XiZa)y2 + fax* - XiZ 2 )y 3 = (7) 

fez/3 - 2/2^3)21 - (xiVs - y&z)zt + (zi?/2 - yix 2 )z 3 = 0. (8) 

In Fig. 4 twice the area (.A 23) of the triangle 023 may be related to the 
coefficient of Xi through the angle (n) between a normal to the orbit 
plane and the x axis. The following equation expresses this and is derived 
in Appendix A : 

2/223 - 222/3 = A 2Z cos n. (9) 

Similar equations may be obtained in like manner for triangles 012 and 
013. When these are substituted into (6), (7) and (8) one obtains 

A23Z1 - A n x 2 + A12S3 = (10) 

A n yi - A u y 2 + A u ya = (11) 

A 23 zi — A 13 z 2 + A X2 z 3 = 0. (12) 

Normalizing to An , (10), (11) and (12) become 

m-iXi — x 2 + 7113X3 = (13) 

W12/1 — 2/2 + W32/3 = (14) 

m\Z\ — zi + W3Z3 = (15) 

where 

wii = A23/A13 (16) 

and 

m z = A12/A13. (17) 

Provided mi and ra 3 are known, we may solve (13), (14) and (15) for 
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Xi , iji and 2, (i = 1, 2, 3). Alternatively we may express Xi , yi and Zi 
in terms of the slant ranges, p, , of the satellite from the ground station 
and solve for these ranges in terms of mi and m 3 . This latter route will 
be pursued. 

Each satellite position measured from the ground station in terms of 
azimuth and elevation may be expressed in terms of right ascension and 
declination, knowing time, by standard techniques 5 (see also Appendix 
B). The direction cosines 15 of the line of sight from station to satellite 
then become 

a, = cos 5, cos on (18) 

hi = cos 5, sin a, (19) 

d = sin Si (20) 

where 

a, , S, = the topocentric right ascension and 
declination of the satellite. 

Since the satellite is essentially orbiting the mass center of the earth, 
it will become advantageous to relate any topocentric coordinates of the 
satellite as measured from the station to the xyz geocentric ones previ- 
ously specified. To do this we write 

Xi = dip, 4- Xi (21) 

y, = bipi + Yi (22) 

Zi = c iP i + Zi (23) 

where 

pi = the slant ranges from the station to the satellite 
Xi , Yi , Zi = the XYZ coordinates of the station in the geo- 
centric coordinate system of Fig. 4. 

For any instant of time the Xi , F, , Z, coordinates of the station may 
be easily derived by referring to Fig. 5. By inspection we write the di- 
rection cosines, which are quite similar to (18), (19) and (20), as 

X, = ft cos Z, cos a,- (0 (24) 

Yi = R cos L sin «,•(*) (25) 

Zi = ft sin L (26) 

where 
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Y L =RCOSLSINO L 




Fig. 5 — Ground station coordinates. 

L = station geocentric latitude (see footnote at end of Sec- 
tion IV) 
R = distance of station from center of the earth (see Ap- 
pendix B) 
at(t) = the right ascension of the station at the instant (t) of 
satellite observation (local sidereal time at the sta- 
tion). 

The local sidereal time* for the above equations may be computed from 

«*(«) = «*»«) ~ Wi (27) 

where 

atig(t) = sidereal time at Greenwich in degrees t 
W\ = west longitude of the station. 

A substitution of (21), (22), and (23) into (13), (14), and (15) will 
accomplish the desired result of expressing .r, , ?/,• , and 2, in terms of the 
slant range and topocentric position of the satellite, the station coordi- 
nates, and mi , m 3 . We obtain 



airriipx — a^pi + dzm^pz = —miXi + X 2 — m 3 X 3 



(28) 



* See Ref. 14, p. 41. 

t The term a, (0 may be computed from the standard meridian time for the 
station by standard techniques. 5 
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bmipi — b 2 p2 + b 3 m 3 p 3 = —miYi + F 2 — m 3 Y 3 (29) 

C1W1P1 — c 2 p2 + c 3 m 3 p3 = —miZi -\- Zz — m 3 Z 3 (30) 

where the unknowns are Wi , m 3 and p, . A solution for p 2 may be obtained 
if mipi and m 3 p 3 are eliminated from these last three equations. This 
yields 

P2 = D/E (31) 

where 

D = otjF, - F 2 + m 3 F 3 (32) 



tf = 



Fj - 



Fi = 



F 3 = 



fll 


02 


a 3 


61 


62 


63 


Ci 


C2 


Ci 


aj 


X, 


(h 


b, 


F, 


h. 


C'l 


Zj 


c 3 


"1 


A' a 


O.i 


bi 


F 2 


&3 


Ci 


Z, 


C3 


«i 


x 3 


a 3 


fa 


r. 


b 3 


t'i 


z. 


c 3 



(33) 



(34) 



(35) 



(30) 



Examining (31) through (36) makes it obvious that p 2 is determined as 
soon as m\ and m 3 are expressed by some separate consideration. 

The one piece of data which has not been worked into the analysis so 
far is the time of flight of the satellite from point 1 to 2 to 3 (see Fig. 4). 
This information added to the above equations, developed solely from 
the geometry of the situation, will permit the solution for the range. 

In Appendix C, it is shown that /Hi and m 3 are not only ratios of tri- 
angular areas formed by the radius vectors to the three satellite positions 
and their respective chords, but also functions of time expressible as 



mi = 



Til 



1 j_ T a( T 2 + ti) , t 3 (t 3 + T1T3 — ti ) dr 
6r 2 3 4r 2 4 dr 

Similar expressions are found in Ref. 9, p. 52, Eq. (5.3). 



(37) 
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A 12 

m 3 = -7- 

" 2 2 1 (38) 

T3 f, ■ Ti(T2 + T 3 ) _ Ti(ti 2 + 7\Tl — T 3 ) rf?' | 

~ r 2 |_ 6r 2 3 4r 2 4 dr "J 

where 

Tl = fc(fc - fe) (39) 

r 2 = *(«, - fa) (40) 

T 3 = Jb(fa - fa) (41)* 

r 2 = geocentric satellite distance at time fa 
ti = time of ith observation 
k = constant for earth satellites f 
= 0.07436574 mm" 1 . 

In (37) and (38) the geocentric range derivatives at this point are un- 
known, and so we approximate mi and m 3 by considering only the first 
two terms in the above equations. We shall call these estimates n Y and 
n 2 respectively. Returning to (31), the range p 2 now may be explicitly 
expressed in terms of known quantities and the unknown r 2 as 

P2 re (l/E)lthFi - F 2 + n 3 F 3 ] (42) 

tt (1/E)[(t V 2 + pu/r^Fi - F 2 + (t 23 + pra/WVs] (42a) 



where 



TV1 



= t\/ti , r 2 s = T3/t 2 (43) 



riT 3 (l + T12) TlT 3 (l + T2») 



TlT 3 \J- T 7-12; m _ '1'3W -T '23/ ( ±A \ 

P12 = g , P23 = g . 144; 

The following geometric relationship between p 2 and r 2 may be written 
by inspecting Fig. 6 

r 2 2 = R 2 - 2R P2 cos tA 2 + P2 2 (45) 

where 

\f/ 2 = the angle at the station between a line to the earth's 
center and one to the satellite 
= 90° + elevation of satellite above a spherical earth of radius 
R. 



* Using here notation similar to Gibbs solution of Ref . 12, p. 100. 

f The constant k is developed in a manner similar to the Gaussian constant for 
the sun, starting with the expression for the satellite period as P — 2rk{a/r,)i, 
where a = the semimajor axis of the orbit ellipse and r e = the equatorial radius of 
the earth. 
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Fig. 6 — The station sphere and subsatellite points. 



Equations (42 ) and (45 ) in combination permit a solution for p 2 and r 2 . 
If one actually proceeds to substitute (42) into (45), an eighth-order 
equation similar to that derived by LaGrange results. This was later 
expressed by Gauss as a trigonometric, transcendental equation solved 
by trial-and-error methods. Rather than pursue this route, it is simpler 
for machine computations to assume an r 2 value (certainly greater than 
R2), calculate a corresponding p 2 using (42), enter (45) with the r 2 -p2 
pair, determine the departure from equality in this equation, increment 
the r 2 estimate, and repeat until the departure from equality is sufficiently 
small. By proper choice of increments, the entire procedure converges 
rapidly to p 2 , r 2 values satisfying the angle-only observations, gravita- 
tional theory, and the approximate mi , m 3 quantities used. The remaining 
p values easily follow from (28 ) , (29 ) and (30 ) . 

In truncating (37 ) and (38 ) after the second term, the range deriva- 
tives were omitted in the above analysis. It is possible to return to these 
equations with these derivatives, which are now crudely known as 
first-estimate averages over the U — U interval. One could then calculate 
second m^ , m 3 estimates from (37) and (38), determine new r< values, 
compute new geocentric range derivatives, and iterate to an n , p, set of 
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solutions. This is roughly the path followed by LaGrangc. By using 
Gauss' relationship involving certain areas of the orbit ellipse, it is pos- 
sible to proceed directly to the true values of ra x and m 3 . This latter 
course of action, which also allows the approximations to be as close as 
needed, will now be followed. 

In Fig. 4, we define as yi the ratio of the areas {Si) formed by the 
radius vectors r x , r 2 , r 3 and the orbit ellipse with the areas {Ai) formed 
by these same vectors and the chords connecting the three satellite 
positions. Thus the sector-triangle ratio 

$i = S 23 /A 23 (46) 

m = Su/Au (47) 

ft = So/An • (48) 

Certainly, in terms of #* , mi and ra 3 may be defined as [see (37) and (38)] 

$232/2 



mi = 



Sny 



(49) 



m 3 = f4 2 . (50) 

If we assume pure Keplerian motion, the radius vector from the earth's 
center to the satellite must always sweep through equal areas in equal 
times. 14 This permits the direct substitution of time differences for the 
S factors in (49) and (50) to yield 



Wx = 



m 3 = 



(fa - t 2 )y 2 = rij/2 ( 51 ) 

{tz — h)yi r 2 yi 

{k — h)y 2 r 3 y 2 



{h — h)y 3 r 2 y 3 



(52) 



It is possible to determine y t knowing the geocentric distance of the 
satellite and the geocentric angle through which it traveled in the meas- 
ured time. This technique will be developed specifically for y x . The re- 
maining variables y 2 and y 3 are derived in identical fashion. Considering 
then the second and third observed positions of the satellite, we write 

« i a k{t 3 - h) Va{l - e 2 ) /-o\ 

r»r 3 sm 2d 

The numerator of this equation expresses Keplerian motion as the radius 
vector sweeps area £23 while the denominator gives twice the area of the 
triangle 023. The variable d is simply the difference in the true anomalies 
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of the satellite at positions 2 and 3. We may shorten the notation by 
writing 

_nVS (54) 

r 2 r 3 sin 2d 

The only unknown in this equation is p, which is a function of the size 
and shape of the final orbit ellipse. Since estimates of r-> and ra arc avail- 
able at this point, p may be eliminated. 4 This permits iji to be expressed 
by a cubic equation* given below: 

Vi - m - aft - a/9 = (55) 



where 



a — 



5/6 + c + e 



(56) 



D v = 2vVsCosd (58) 



c = 



r 2 + n - D y 



2D„ 



(59) 



52 3 



'-i5*' + ws x + - (00) 

x = b/yi 2 — c. (61) 

It is indeed unfortunate that the factor a is ultimately a function of y x , 
for if this were not the case an explicit solution for iji would exist. It 
turns out that e is sufficiently small that it may be equated to until a 
first estimate of y\ is found, whereupon it may then be calculated from 
(61) and (60) and used to derive a better estimate of yi . The estimate 
of iji may be improved as much as desired by traversing this loop a 
multiplicity of times. Knowing the y t values, w?i and m 3 are next com- 
puted from (51) and (52) and entered in place of the approximate n\ 
and n 3 values in (42). The three ranges to the satellite are then cal- 
culated as described earlier. 

VI. GENERATING THE ORBITAL ELEMENTS 

6.1 General Procedure 

Three angular positions of a near-earth satellite as observed from a 
given ground station are sufficient to generate a set of orbital elements. 

* See also Ref . 9, p. 56. 
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The addition of range data permits elements to be computed which arc 
generally more accurate than those developed from angle only informa- 
tion. It will now be shown how these elements may be generated first 
without regard to any perturbations — as if the satellite orbited a spher- 
ical, homogeneous earth. The perturbations arising out of the earth's 
oblateness will then be computed, and finally their effects will be appro- 
priately combined with the input data to produce the modified orbital 
elements described in Section III. 

6.2 Generating the Unperturbed Elements 

6.2.1 Locating the Subsatellite Point. 

Given the azimuth, elevation, and range of a satellite from a specified 
ground station, we begin by calculating the latitude-longitude coordinates 
of the so-called subsatellite point. This is the point of intersection of a 
line, drawn from the satellite to the earth's center, with the surface of 
what we shall call the station sphere. This sphere is centered on the 
earth's center and has a radius that equals the geocentric distance of the 
ground station. Fig. 6 shows this. Here the central angle g% , which lies 
between the station A and the subsatellite point *S', is computed by 
solving triangle AOS. To do this, we recognize the angle OAS is known 
from the measured elevation of the satellite reduced to a spherical earth 
(see Appendix B). Thus 

ZOAS = 90° + E (62) 

where 

E = elevation in degrees of the satellite above a plane tangent 
to the station sphere at the station. 

Since the slant range to the satellite, p, and the geocentric distance of the 
station (see Appendix B) are also known, we write 

cos (E + g 2 ) = (R/p) sin g 2 (63) 

from the law of sines applied to triangle AOS. Solving for g 2 we obtain 

■-H»nH*J- (64> 

Since g 2 must always be less than 90°, there is no angular ambiguity in 
this quantity as expressed above. 

From the same triangle, the geocentric distance of the satellite becomes 

r = pCOsg . (65) 

sin 02 
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The latitude and longitude coordinates of the subsatellite point are 
obtained by solving the spherical triangle ANS knowing g 2 and the azi- 
muth of the satellite, (A). By the law of cosines, the latitude of S be- 
comes 

L = sin -1 [cos g 2 sin L a + sin g 2 cos L a cos A]. (66) 

Again, there can be no angular ambiguity since L ranges from —90° to 
+90°. 

The longitude of S' is determined from the longitude offset (6) of 
S' from A . By the law of sines 

6 = sin- r sing2 7 L l . (67) 

|_ sm A J 

Angle b ranges from —180° to 0° to +180°, but the direction (sign) of 
this longitude offset is directly determinable from the azimuth of the 
satellite as seen from station A. Thus the subsatellite point has been 
located in latitude and longitude on the station sphere.* These corre- 
spond to the geocentric latitude and longitude of »S' measured on the 
earth. 

6.2.2 Computing Orbit Inclination and Ascending Node 

Two subsatellite points must be known for the geometrical calculation 
of the orbit inclination and the ascending node. Consider Fig. 7, where 
Si and S 2 represent two known positions of the subsatellite point. The 
orbit inclination angle i is calculated in the following manner. The great 
circle arc g i2 is first computed by standard techniques, knowing the lati- 
tude and longitude of Si and S 2 . It is quite important, however, to re- 
member that the earth has spun through a specific angle during the time 
that the subsatellite point progressed from Si to S 2 , and both nodal 
regression and apsidal advance have occurred. Therefore, S 2 must be 
assigned a new longitude suitable to the nonrotating geometry of Fig. 
7. This new longitude in degrees is 

I, = U - %^ 360 + N (68) 

where 

l 2 = west longitude of S 2 from the Greenwich Meridian 
computed from the satellite azimuth, elevation, 
and range 

* Certain special rases exist, such as when S' lies on the same meridian as 4. 
These require no additional data and are actually simpler, since then b = 0, g< = 
I j — A,,, and spherical triangle ANS' is not needed. 
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]f jg, 7 — Geometry for orbit inclination and ascending node calculation. 

I,, — h = time in minutes between the occurrence of the two 
measured subsatellite points 
1436 = time for one earth rotation in minutes of mean solar 
time 
N = nodal regression in time U. — h • 

We shall not compute the effects of apsidal advance upon the lati- 
tude and longitude of the subsatellite points, as this perturbation will be 
taken into account in the following section by adjusting the true anoma- 
lies. 

From Fig. 7, one writes 



sin i = 



sin Li 



sin Li 



sin c sin (c + <7i 2 ) ' 
Using the last two terms of (69), one may solve for c to obtain 

sin #12 



c = tan 



sin L 2 
sin h\ 



— COS 012 



(69) 



(70) 



Here sin gn is always positive, since Si and So correspond to satellite 

* N is initially set equal to zero until an unperturbed set of elements is gener- 
ated from which A T may be determined: see Section VII. 
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positions in sight of station A. The arc c, however, may range from 0° 
to 180°. If the denominator of (70) is positive, c is between 0° and 90°; 
if negative, c is between 90° and 180°. 
Knowing c, i is simply 



[sin Li~| 
sin c J 



t-HiT 1 =£=M (71) 

|_ sin c J 

with no ambiguity, since it ranges from 0° to 90° for proper satellites. 
In similar fashion, b of Fig. 7 is 



b = cos 



<n * C -]. (71a) 



cos L 



If cos b is positive, b is between 0° and 90°; if negative, between 90° and 
180°. The ascending node is then simply 

fii = /i - b ± 360 n (71b) 

where 

h = the W longitude of Si. 

n = an integer chosen to cause S2i to lie between 0° and 360°. 

So far only the geometry for »Si and S 2 both on the same side of the 
equator has been considered. If they are on opposite sides, then the only 
change is in (70), where the — sign becomes a + sign. 

6.2.3 Computing the True Anomaly along with Orbit Size and Shape. 

If the latitude and longitude positions of the three subsatellite points 
are known along with the corresponding geocentric satellite distances, 
the orbit size and shape may be computed and the satellite's true anomaly 
determined. To do this, consider Fig. 8, where n , r 2 , r 3 , v 2 / and v u are 
all known. The values for n , r 2 , r 3 were obtained by (65), and v n , v n 
are merely the great circle arcs between S 2 and Si and between S3 and 
Si, respectively, corrected for earth spin and later for regression per 
(68). It is assumed that perigee moves in a secular manner along the 
orbit ellipse during the time of travel of the satellite from 1 to 2 to 3. 
Consequently the true anomaly differences determined geometrically 
must be corrected as follows 

Vn = «>2i + ^ (72) 



2 



7T 



v 3 i = v tl + ^ (72a) 
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FOCUS OF ELLIPSE 
(EARTH CENTER)" 



\ 


V 31 


\ \ 



v' 



Fig. 8 — True anomaly geometry. 



where 

^21 , ^23 = the geometric true anomaly differences 

AP = the apsidal advance in degrees per anomalistic 
period 
Mi , Mi , M A = the mean anomalies* of the satellite at 1, 
2, and 3 respectively. 

Quite clearly the mean anomalies are not known until the position of 
perigee is determined. Thus, the corrections of (72) and (72a) are not 
applied to y 2 i and v 3 i until first estimates of the perigee position and of 
the apsidal advance rate are obtained. At that time a first estimate of 
the true anomalies of the satellite is known, and this is used to deter- 
mine a first estimate of the mean anomalies. With these values, (72) 
and (72a) are used to correct the purely geometrical true anomalies and 
the loop is traversed until no further significant changes in v 2 i or y 3 i are 
observed. 

Let us now set 



* See Ref. 9, p. 113. 
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p = a (I - e) (73) 



where 



a = the semimajor axis of the orbit ellipse 
e = eccentricity of the orbit ellipse. 

Then for the three satellite positions the polar form of the ellipse* be- 
comes 

n = n^ (74) 

1 -+- e cos ih 

r, = — - V , - . (75) 

I -\- e cos (vi + v-n) 

r » = T-T V ( ■ x CW) 

1 + e cos (wi + y 3 i) 

whore 

i>i = the true anomaly of the satellite at position 1. 

It is helpful to rewrite these three equations as 

e cos Vi = (p/ri ) — 1 (77 ) 

1 + Ac cos vi + Be sin Vi = p/r? (78) 

1 + Ce cos vi + De sin ^i = p/r 3 (79) 

where A , B, C and /) are all known quantities defined as 

A = cos v n (80) 

B = -anwa (81) 

C = cosy 31 (82) 

D = -sin v 3 i . (83) 

Solving for p in (77), (78) and (79) yields 

_ C - 1 + DQ. - A)/B 
P "C_1,D/1 A\' (84; 

T\ r 3 



1 + 2(1 _4V 



All quantities in the right-hand member of this equation are known, so 
that p is determined. Returning to (77) and (78), the solution for the 
true anomaly is 

* See Ref. 9, p. 110. 
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-i [" sin fi ~1 
|_COS ViJ 



Vi = tan * I '^T-l± I = tan 1 



B r 2 



Ut-i 









(85) 



unless (p/ri) — 1 is zero, which occurs for a circular orbit. In this case, 
the solution for v x is of little importance and may be arbitrarily set to 
90°. 

In general, the true anomaly ranges from to 360°, but all quadrant 
ambiguities in (85) may be resolved by observing the algebraic sign of 
sin Vi and cos v t , expressed by the numerator and denominator of the 
final fraction in (85). 

With vi known, the eccentricity from (77) is simply 



t - 

rx * (86) 

e = . 



£ - 1 



COS Vi 

The radius of perigee from (73 ) is 

*-'(£$)■ (87) 

6.2.4 Reduction to Epoch 

One of the features of the modified orbital elements is that they are 
quoted for the epoch at which the satellite passes through perigee. The 
node and perigee positions as well as the true anomaly of the satellite 
determined in the preceding sections were computed for the time of the 
first data point. They must now be reduced to the time of perigee passage. 

During the time the satellite moved from perigee to the first data 
point, the earth rotated, the orbit node regressed, and the perigee moved. 
The angular velocity of the earth is known and the node and perigee 
motions are secular functions of i, e, and r p , as will be shown later in 
Section VII. Consequently it is possible by calculation to "move the 
satellite" back to perigee, taking into account these other motions as 
well. To do this we shall assume Keplerian motion along an ellipse whose 
perigee moves in time. If we knew the true anomaly of the satellite at 
position 1, referenced to the moving perigee, we could determine the 
time required for the satellite to orbit from perigee to data point number 
1. This true anomaly may be obtained in an iterative manner using 

vi = vi + AP(Mi/2t) (88) 
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where 



AP = apsidal advance in degrees per anomalistic period 
Mi = the mean anomaly of the satellite at data point 1 
V\ = the geometrical true anomaly of satellite when at point 
1 (from Section 6.2.3). 

A first estimate of M i is obtained from Vi using Kepler's equation. This 
permits a first estimate of V\ by using (88), which in turn allows a sec- 
ond estimate of M\ based on the first vi estimate. The procedure cycles 
until further changes in V\ are insignificant. The final value of M\ corre- 
sponding to the final v x is used to determine the time of flight of the 
satellite from perigee to the first data point, which is 

h = {Mi/2*)T a (89) 

where 

T„ = the anomalistic period in minutes. 

The elements arc back dated by the amount h to produce modified orbital 
elements. In these the ascending node is simply 

& - « x + AQt b (90) 

where 

S2i = the west longitude of the node at data point no. 1 
Aft = the regressional rate in degrees per minute. 

Similarly, the argument of perigee is 

P = Px+ &P{Mi/2tt) (91) 

where 

Pi = the perigee position at data point no. 1 
AP = the apsidal advance rate in degrees per anomalistic 
period. 

6.2.5 The Anomalistic Period 

The anomalistic period may be computed from the mean anomalies 
at the data points (after these are known from the first estimate of the 
orbital elements) and the times of flight through these points. For ex- 
ample, one may write 

2tt(< 2 - h) , . 

la Mi -Mi ^ 2) 
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or 



or 



M 3 - Mi 



M 3 — Mi 

Equally well, T a may be computed from the finally determined values 
of i, e, and r p by methods outlined in Section VII. 

6.3 Adding the Perturbations to the Orbital Elements 

In the derivation of the modified orbital elements, perturbations seem 
required before sufficient orbital data are available for their computa- 
tion. The general technique is to assume all perturbations zero until a 
first set of elements has been calculated for epoch t t , then to supply 
these perturbations to produce a second set of elements, and to continue 
this iteration until further element changes are insignificant. Thus the 
N of (68), the AP of (72), (72a), and (88), and the T a of (89) follow 
this route. 

For orbits of the TELSTAR satellites which have had eccentricities 
of up to 0.4 and r p values as low as 4556 statute miles, convergence by 
the above method occurs within three iterations. Four iterations are 
normally required for eccentricities up to about 0.8. 

VII. THE SECULAR PERTURBATIONS 

7.1 Perturbation Theory and Formulae 

The perturbations to be outlined here are secular and produced only 
by the earth's oblateness. The gravitational field therefore will be inde- 
pendent of the longitude and symmetrical about the earth's equator. Its 
potential* may be expressed as 

GM [l - t Jn (^)" Pn(C0B L c )] (95) 



u = 

r 



where 



J n = constants determined in the first instance by ob- 
servation of many satellites 



* See Ref. 16, p. 4. 
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R E = the equatorial radius of the earth = 3963.347 statute mi. 
r = the geocentric satellite distance 
P„(cos L e ) = the Legendre polynomial of order n 

L c = the geocentric colatitude of the satellite. 

Writing (95) using only the J 2 terms, we have 

U = ?M. [l + I J 2 (^J (1-3 cos* L c )] (96)* 

where 

GM = Newton's universal gravitational constant times the 
mass of the earth 16,17 
= (3.44280 ± 0.000026) X 10 8 statute miles'/minutes 2 
J 2 = f x Jeffreys' ./ (see Refs. 18 and 19) 
= 1.0827 X 10~ 3 . 

Equations of motion for a satellite orbiting this oblate earth are found 
by setting the spherical polar components of acceleration equal to 

du idU . i du 

a" I " KT> aild ~ r-a-i 

dr r dL c r sm L r . d<p 

where r, L c , and <p are the spherical coordinates of the satellite. After a 
bit of work, the following secular rates of the orbital elements are pro- 
duced. For the secular nodal regression f in radians per nodal period we 
have 

A£2 = -3ttJ 2 ( — Y cos i + OJ 2 



J «(t)' 

(97) 



= -0.16028 ^-'lO 6 
P 2 



where 



* J i = since the orbit plane is assumed to include the center of mass of the 
earth. 

t Many authors have derived perturbation equations for orbital elements and 
expressed them in numerous ways, quite often with great elegance and detail. It 
may seem, therefore, that (97) is unduly simple. This analysis, however, is con- 
fined to secular rates only. It can be shown that contributions to these rates from 
the third through the sixth zonal harmonics of the earth's potential function (as 
well as from the J 2 2 terms) are about three orders of magnitude less than the J« 
contribution given by (97). They are, therefore, omitted in this analysis where 
long-term predictions of pointing angles need only be accurate to within 0.05°. 
Secular contributions to the orbit eccentricity and inclination are also of order /2 s 
and omitted. (See Ref . 21 for expressions of higher-order terms.) 
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p = a(l — e ) 
OJ 2 = terms of the order of J 2 which are omitted. 

The secular apsidal advance rate becomes 

Ac* - - 1 7rJ 2 (— Y [1-5 cos 2 * + Oe 2 ] 

= 0.08014 (5cos^--l) ]()6 
P 2 



(98) 



in radians per nodal period. Finally, the nodal period which expresses 
the energy in the ellipse is 

^A)T i -i^(f) ,(g ^^ i)(i -^] (99) 

T 2 J o ! (l - e 2 ) 3 

where 

2v/(GM) { = 3.38629354 X 10~ 4 minute/statute miles 3 ' 2 . 

The prime sweep interval of the modified orbital elements is related 
to the true nodal regression by the following expressions 

Aa , _ AO X 82505.92 deg/day (9fc) 

■I n 

PSI - 1440 - 4(AQ/ + °- 9856) minute, (qqh) 

Atf + 0.9856 ("W 

+ 360 

The rates for the modified orbital elements are not complete until the 
anomalistic period has been derived. This may be accomplished know- 
ing T n and Ato. The procedure is outlined below. 

7.2 Determining the Anomalistic Period from the Nodal Period 

We have denned the anomalistic period as the time required for the 
satellite to move from perigee to the next perigee. The two perigee posi- 
tions are shown as points pi and p 2 in Fig. 9, where it is seen that the 
angular motion of the satellite in the time T a may exceed 360°. During 
this time the perigee itself is in motion and its rate is expressed by (98). 
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ANGULAR MOTION INDICATED 

TAKES PLACE IN TIMES 

T n OR T a AS SHOWN 



Fig. 9 — Relationship between T n and T„ . 

Thus the anomalistic period exceeds* the nodal period by an amount t 
so that 



or 



T a = T n + t 

r a = (i +/)r„ 



(100) 

(101) 

(102) 



where 



t = the time for the satellite to move from pi to p 2 
/ = the fraction of T» represented by t 
M = the mean anomaly corresponding to /. 

We may also express the motion of perigee in time T a as the angle 

v p = Aco TJT a . (103) 

If T a is initially set equal to T„ , then (103) represents a first estimate of 
v,, from which a first estimate of/ may be computed from 

* If Aw is negative, T n exceeds T a , but the conversion equations are still valid. 
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£ = 2tan -.[(L^i) , tan |] (104) 

and 

M = E - esiaE. (105) 

By (102) a first anomalistic period estimate is obtained which can yield 
a second estimate of v p by (103), and the process is iterated until succes- 
sive changes in T a are sufficiently small. Since (102) through (105) con- 
tain essentially two unknowns, v p and T a , closed-form approximate 
solutions are possible by using truncated series of expansions of the 
tangent and sine functions. The iterative method, however, has the 
advantage that the accuracy is limited only by the number of times that 
the user is willing to traverse the loop. 

Knowing the anomalistic period, it is straightforward to express the 
nodal regression and apsidal advance in degrees per anomalistic period 
suitable for use in the modified orbital elements. 

VIII. ROUND-THE-WORLD ORBIT PREDICTIONS 

So far, techniques have been presented for developing orbital elements 
suitable for predicting future satellite pointing angles. These predictions 
suffer little error when referred to the ground station which took the 
original data. However, if predictions are made for other stations, the 
errors can be greater due simply to the fact that the elements were de- 
rived from only a small portion of the orbit ellipse — a fact which makes 
it difficult to assess the size and shape of the ellipse. The effect will be 
most pronounced if the latitude of the second station differs greatly 
from that of the first* as, for example, the case of elements generated 
from Andover, Maine, trajectories and used for Woomera, Australia, 
predictions. In previous works this problem is simply solved by obtain- 
ing data from stations spaced around the world and using it in orbit 
generation. In this paper it is desired to confine the data taken to one 
station and yet predict for remote stations located anywhere. This is ac- 
complished by recognizing that, during the course of a day, a single base 
station sees a good portion of the orbital ellipse on successive passes. 
For each such pass a set of elements may be generated. Each set will 
about equally well describe the orbit over the base station, but can di- 
verge in predictions for stations at different latitudes as illustrated by 

* It should be recognized that stations at different longitudes (but at equal 
latitudes) are not to be considered unique, since they essentially observe the same 
portions of the orbit at different times of day. 
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exaggeration in Fig. 10. This is due principally to estimating the size 
and shape of the orbit ellipse from data which necessarily contain some 
measurement error. One could improve the situation at remote stations 
by placing an ellipse of best fit through all the essentially different base 
station trajectories and thereby gain a better estimate of the eccentricity 
and radius of perigee in the orbit. If base station predictions made from 
elements that fit each individual pass are excellent, one could equally 
well fit an ellipse by a least-squares method to the predicted rather than 
the measured base station data, and still improve the situation at remote 
stations. As a matter of fact, proceeding from good base station predic- 
tions tends to smooth the input data by bringing them into accord with 
the motions permitted by gravitational theory and this, in a sense, filters 
some of the observation noise. 

The technique for this round-the-world fit proceeds in detail as follows. 
Modified orbital elements are used to predict the true anomaly and the 
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corresponding geocentric distance of the satellite for a set of passes seen 
from the base station. The polar equation of the ellipse is then written 
as follows : 

1 1 e cos v finR\ 

r a(l — e 2 ) a{\ — e 2 ) 

= B + Bxcosv. (107) 

The predicted true anomaly-geocentric distance pairs are then used to 
determine B and Bi of (107) by standard least-squares methods. It is 
important to note that the true anomaly v of (107) must always be re- 
ferred to the moving orbit perigee as was done in Section 6.3. 

Providing that e and r p do not seriously change then magnitude during 
the set of passes being considered (that is, no great orbit energy loss), 
we may determine these quantities from the least-squares B and Bi 
values to be 

e = Bi/Bo (108) 

"-sn^V (109) 

While the determination of orbit size and shape for round-the-world 
orbit predictions may be carried out by least-squares methods, the 
node position and argument of perigee are often treated differently, since 
they are not determined to the same degree of accuracy by each pass of 
the satellite over the base station. Pass selection is often a good procedure 
to follow for these parameters. For example, data from a pass that con- 
tains perigee are generally considered much more likely to produce a 
more nearly accurate measure of the perigee argument than those from 
a pass which does not. The same is true of the ascending node. 

IX. IMPROVING THE ELEMENT RATES BY DIRECT MEASURE 

Quite obviously, if the modified orbital elements are known at two 
separate epochs, the secular rates may be computed directly. Let us 
assume that the elements are correct except for the rates of nodal regres- 
sion, anomalistic period, and apsidal advance. 

9.1 Nodal Regression 

The west longitude of the ascending node at the earlier epoch (fii) is 
related to its position at the later epoch (fi 2 ) by 

ih + 360A - ih = ^^ degrees (110) 

rbl 
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where 

t n = time interval in minutes between the two epochs 
PSI = the nodal regression expressed as the prime sweep in- 
terval (1440 — M minutes) 
A = an integer which may have any positive, negative or 
zero value. 

Examination of (110) reveals ambiguities in multiples of 360° in the 
determination of the prime sweep interval, since without some knowledge 
of its magnitude, one cannot be sure if the node has regressed more than 
one revolution to its fl 2 position. The theoretical value of the prime sweep 
interval is available, however, by the techniques of Section VII. Thus the 
ambiguity is avoided unless t 2 is extremely large, which allows even small 
variations in PSI to promote as much as one revolution difference in 
nodal regression. Solving (110) for the prune sweep interval, we have 

PSI - o o^zW minutes - < m ) 

9.2 Anomalistic Period 

From the epochs of the two sets of modified orbital elements on hand, 
a single anomalistic period, having zero rate of change, may be com- 
puted. First, the whole number of satellite passages through perigee in 
the time t n between epochs is expressed as 

A n = INTF {t n /T at ). (112) 

Then a refined estimate of the anomalistic period becomes 

T a = tn/An (113) 

where 

T at = theoretical anomalistic period computed from the first 
set of MOES using (99) through (105) 
INTF = indicates that only the integer portion of the division 
is to be retained. 

Sometimes the anomalistic period of a particular set of elements has 
been well refined by the above procedure, so that it is quite precise over 
a short period of time. If, some time later, a second set of elements is 
determined and shown to fit the orbit well at this later time, then it is 
often a good procedure to retain the initial anomalistic period and com- 
pute a period change to be consistent with the satellite's arrival at perigee 
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on epoch 2. This incidentally, is the only manner in which a change in 
period is deduced for the TELSTAR satellites. The procedure is as fol- 
lows. Calculate A„ using (112) as before. A first estimate of the change in 
period per period becomes 

ATa = 2tn ~ j nTat . (114) 

A n 2 

It is possible that the A„ value calculated by (112) is not correct owing 
to the change in the period during the interval tn . Therefore we compute 
a time interval tf as a trial to compare to t i2 

t f = A n T at + (A„ 2 /2)AT a (115) 

and if t f exceeds t vl , we decrease A n by 1 and repeat the operations of 
(114) and (115) in an iterative sense until t f is less than t n . (Of course, 
if tf ever equals t a in the process, then the corresponding A n value is cor- 
rect and no further iteration is needed.) When t f is less than t u compute 

t n = (An + l)T ot + |(A« + 1) 2 A7' U (110) 

and if t n exceeds t i2 , the current A n value is valid. If tn is equal to t n 
then A„ should be decreased by one. If t n is less than t n , increase A,, 
by one and repeat the process from (114) on. Eventually the procedure 
converges on an A„ value that either makes tf equal to or less than t vi 
and tn greater than t vl . The corresponding A2'„ is the proper one to use. 

9.3 Apsidal Advance Rate 
The apsidal advance rate may be expressed as 

Aco = "* " ** + 360B (117) 

An +f 

where 

, m fa - t f ( U8 ) 

J Tat 

and 

B = an integer chosen by comparing Aw with the theoretical value 
computed from the first set of elements, as was done with A of 
(110). 

In cases where no AT a is computed, the A„ value of (112) and / = 
should be used in (117). 

The advance rate of (117) is expressed in degrees per anomalistic 
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period, which has direct meaning provided no change in period is cal- 
culated. Whenever AT a is not 0, (117) produces the advance rate which 
when applied to each new period* between epoch 1 and epoch 2 results 
in the motion of perigee from a>i to co 2 as stated in the elements. 

X. A COMPUTER-OPERATOR ENSEMBLE FOR ORBIT GENERATION 

Computational techniques have been developed to permit modified 
orbital elements to be determined initially from the powered flight 
parameters of the launch vehicle and subsequently from measured satel- 
lite trajectories over a ground station using techniques previously de- 
scribed in this paper. These have been augmented to enable the elements 
to be refined so as to make them valid for several months and to monitor 
their validity on a routine basis. The techniques have been incorporated 
into five computer programs which have been routinely used, in the 
case of the second TELSTAR satellite, to produce elements for antenna 
drive tapes, pass scheduling, and satellite attitude predictions. It is the 
purpose of this section to show the means by which the programs are 
integrated into a computer-operator system. 

10.1 The Computer Programs 

The five computer programs employed for orbit generation, refine- 
ment, and monitoring techniques are: 
(t) The MOEGEN G3 program 
(it) The MOERATE program 
(tit) The REFINE PERIOD program 
(w) The CORD DOPPLER G2 program 
(v) The MONITOR TAPE ANALYSIS program. 
The first four of these were written by the author, while the last 22 was 
basically developed by D. A. Aaronson with certain subroutines by D. A. 
Ramos. In this section we shall discuss the general function of each pro- 
gram. 

10.1.1 The MOEGEN G8 Program 

The MOEGEN G3 Program generates modified orbital elements for a 
near-earth satellite given a list of azimuth, elevation, and, if available, 
range of the satellite as seen from a ground station at specified times. If 
just these data are used, the resulting MOES are said to be free-fit. If the 
PSI, apsidal advance, and anomalistic period are known from past data, 
the MOEGEN program may be set to accept these rates along with 
trajectory information as before and produce a suitable set of MOES. 

* A new period is calculated each time the satellite passes through perigee. 
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These MOES are said to be forced. Since the anomalistic period cannot 
be accurately determined from a single pass, this latter mode of opera- 
tion is most valuable in producing elements of long lasting validity. 

In addition to generating elements, the MOEGEN program compares 
trajectories generated from each set of MOES it produces with the points 
used in the generation, as well as with any additional submitted points 
from any of the passes from that ground station. Here the program re- 
ports azimuth differences, elevation differences, and great circle arc 
pointing errors as well as range errors. Since the program produces a set 
of MOES for every group of three azimuth-elevation-range points, these 
comparisons aid in the selection of the best set of MOES from the group 
generated. The actual selection, which will be discussed later, is done by 
an operator, since typically only a few sets of MOES are produced on 
each generation. 

10.1.2 The MOERATE Program 

The MOERATE program, as its name suggests, refines the PSI, 
apsidal advance, and period rates along the lines described in Section 
IX by examining two sets of elements, usually the ones previously in use 
and the one just generated. In addition, it compares the rates so estab- 
lished to theoretical rates determined from the previous set of MOES us- 
ing the gravitational theory given in Section VII. In generating rates, 
higher-order changes beyond the PSI, apsidal advance, and change in 
period are neglected. For the TELSTAR satellites this omission has 
resulted in no serious errors in the orbit. 

10.1.3 The REFINE PERIOD Program 

Because computed satellite trajectories are most sensitive to the value 
calculated for the anomalistic period, the REFINE PERIOD program 
offers an alternate means for its determination and is particularly valu- 
able before monitor tapes* are available. The program requires a set of 
MOES and times of meridian crossings from any specified ground sta- 
tion. From each of these it computes an anomalistic period assuming all 
other element values to be accurate. If there are no errors in the recorded 
time for each meridian crossing, in the adopted station coordinates, or 
in the elements, quite naturally all the computed period values should be 
identical. Any scatter in the computed periods points up errors in the 
above quantitites. If the meridian crossing times and station locations 

* These are magnetic tapes generated at the ground station. They contain, 
among other data, the measured trajectory of the satellite. 
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are known to be correct, the scatter acts as a figure of merit for the 
generated modified orbital elements. 

10.1.4 The CORD DOPPLER G2 Program 

This program accepts either MOES or, in conjunction with a sub- 
routine called MOE, expected burnout parameters for the final powered 
stage of the satellite launching vehicle. From these it produces an ephem- 
eris for passes over any specified ground station. If burnout quantities 
are used, the program converts these to MOES, which it prints prior to 
the ephemeris listing. In addition, this program also produces world maps 
of the satellite ground trace, drive tapes for both the Andover and 
Holmdel ground stations, and reports the solar illumination status of the 
satellite along with a number of communication quantities. 

10.1.5 The MONITOR TAPE ANALYSIS Program 

The monitor tape analysis program written by D. A. Aaronson and 
described in detail elsewhere 22 compares a satellite trajectory as recorded 
on a compressed or uncompressed monitor tape 22 to that produced by a 
set of MOES. Its output is a printed sheet giving great circle arc pointing 
errors and slant range errors (provided ranging data are present on the 
monitor tape). It also indicates received microwave carrier power so that 
one can ascertain whether or not the Andover horn antenna was indeed 
pointed toward the satellite. This is, of course, essential to proper assess- 
ment of the printed errors. 

1 0.2 Orbital Element Generation Operational Procedures 

The generation of modified orbital elements for the TELSTAR satel- 
lites falls into two categories: the initial generation just after launch and 
the subsequent generation after monitor tapes have become available. 
Figs. 11 and 12 show the logic operations for these two cases. 

10.2.1 Initial Element Generation 

The initial elements are determined from the expected burnout pa- 
rameters* of the final powered stage of the launch vehicle. These 
parameters are based upon the telemetry from the first two stages. Basic 



* For the TELSTAR satellites, these have been supplied by J. W. Timko of 
Bell Laboratories. 
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to these parameters is the assumption that the third stage performs in a 
design center manner. 

Subroutine MOE of the CORD DOPPLER G2 program accepts the 
burnout parameters as shown in Fig. 11 and calculates corresponding 
MOES, which are passed to the main program for ephermeris generation 
listings. 

Modified orbital elements calculated from the burnout parameters are 
submitted along with subsequent meridian crossings to the REFINE 
PERIOD program to produce an anomalistic period list as described in 
Section 10.1.3 (see also Fig. 11). The operator examines this list for 
consistency. Because the REFINE PERIOD program bases its period 
calculation on what it assumes as an otherwise correct set of MOES and 
because the final stage of the launch vehicle, often without telemetry 
provision, may have performed in a manner different than assumed for 
orbit generation, it is usually satisfactory to call the period list con- 
sistent if the spread of values does not exceed a few tenths of a minute 
or about 0.1 per cent of the mean anomalistic period for the distribution. 
If the list is satisfactory according to this criterion, the operator substi- 
tutes the mean value of the period into the set of "burnout" MOES and 
checks minute by minute the predicted trajectory with measured ones by 
means of the MONITOR TAPE ANALYSIS program as soon as the 
monitor tape becomes available (see Fig. 11). 

Great circle arc errors and range errors reported by the monitor pro- 
gram are scanned by an operator, and if these are within established 
limits, the element generation process stops and the set of MOES on 
hand is used for future predictions. If limits are exceeded, the operation 
proceeds as outlined by Fig. 12, to be discussed later. 

Care should be exercised in establishing error criteria for this operator 
scanning. The first set of MOES, being based on expected burnout 
parameters, cannot generally produce pointing errors below 0.1° in arc 
or ±10 miles in range. 

If the anomalistic period list originally generated by the REFINE 
PERIOD program has too great a spread, an inconsistent set of MOES is 
generally indicated (see Section 10.1.3). This can occur when the final 
stage of the launch vehicle departs sufficiently from its design trajectory. 
It is useful at this time to have on hand a number of MOES generated 
from burnout parameters with tolerances within the final-stage per- 
formance limits. These can be used with the REFINE PERIOD 
program with the expectation that at least one set will produce an ac- 
ceptable anomalistic period list. If this procedure is unsuccessful, there 
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is no recourse but to wait for the monitor tape or sufficient data from 
tracking stations to establish a correct orbit prediction. 

10.2.2 Subsequent Element Generation 

As soon as monitor tapes become available from Andover, the opera- 
tion passes into what is termed the subsequent orbital element generation 
stage. This is outlined in Fig. 12. 

As before, monitor tape trajectories are compared to those generated 
by the current orbital elements. If these comparisons are unsatisfactory 
(see Section 10.2.1), the horn azimuth, elevation and the measured range 
of the satellite along with the corresponding times of measurement are 
fed to the MOEGEN G3 program. Only a free fit (see Section 10.1.1) is 
made if no previous rate data are available. Both free and forced fits are 
made if reasonably reliable values for the prime sweep interval, apsidal 
advance, and anomalistic period are on hand. In the early phases of 
orbit determination, the REFINE PERIOD program is used to establish 
an anomalistic period by meridian crossings. After the first workable set 
of MOES is established by the subsequent generation procedure, the 
MOERATE program then supplies most of the period values. 

In normal operation, the MOEGEN G3 program will compare the 
input data points to trajectories generated from all output MOES. 
Optionally, past trajectory data may also be submitted to the program 
for similar comparisons. An operator scans these comparisons, particu- 
larly those concerning the trajectory from which the MOES were de- 
rived. The free -fit comparison should have great circle arc differences 
below 0.05° and range errors below a few miles provided the input tra- 
jectory data are consistent with the pointing angles resolved to within 
0.02°, as is generally true on the Andover monitor tapes. The forced 
fit may produce greater errors, which naturally depend on the correct- 
ness of the forced data. By these comparisons the operator selects a 
single set of MOES. Failure to select can occur if input data are ex- 
cessively noisy or for some reason inconsistent. There is no alternative 
in this case but to wait for the next monitor tape. 

The MOES selected are used directly for future predictions if no 
previous set other than those from burnout exists. The burnout MOES, 
being only expected values and having only an approximate epoch, are 
generally not suitable for submission to the MOERATE program. If a 
previous set of elements from MOEGEN G3 are available, they are com- 
pared to the current set by MOERATE to determine more accurately 
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the PSI, apsidal advance, and anomalistic period. The operator sub- 
stitutes these for the earlier rates and again enters MOEGEN G3 using 
the forced fit mode only. The resulting elements which best fit the pres- 
ent and past trajectories are selected for subsequent use. 

XI. TYPICAL EXPERIMENTAL RESULTS 

In this section, a typical orbit generation procedure will be illustrated 
along the lines described in Section 10.2.2 and then a theoretical round- 
the-world orbit fit will be shown using the techniques of Section VIII. 
Data for the typical element generation will come from 20 minutes of 
measurement on each of two passes occurring over Andover on June 30 
and July 30. Data from passes on other occasions will be used to determine 
the accuracy of the predicted trajectories but, for purposes of this il- 
lustration, in no way will these latter data be used to refine the orbital 
elements except as specifically stated. 

Fig. 13 shows an output from the MOEGEN G3 program operated in 
the free-fit mode. Three sets of modified orbital elements were calculated 
from three data points (azimuth, elevation, range) taken on June 30, 
1964. These three sets, shown on lines 7, 8, and 9, differ from each other 
only in the method of determining the anomalistic period. In the first set, 
the period is determined from the elements using the perturbation for- 
mulations of Section VII. In the second and third sets, the period is de- 
termined from the time of flight of the satellite between the data points 
per Section 6.2.5. 

Directly beneath the elements, the program states that the orbit offset 
is 31.48778°. This is merely an indicator which means that the station 
is that number of degrees of latitude removed from a point of equal 
longitude on the orbit plane at the time of the first data point. If this 
quantity were zero, then the station would be in the orbit plane and 
this three-point analysis would be invalid.* 

The portion of the output enclosed in parentheses displays the com- 
parison of predicted -to-measured trajectories. There are three sets of 
comparisons, corresponding to the three separate sets of modified orbital 
elements generated. Here we see that arc and range errors are reasonably 
small for the June 30 pass (which data were used to determine the 
elements) but that the fit on June 10, for example, exhibits 5° pointing 
errors. A survey of the fits clearly indicates that the orbital rates (i.e., 
prime sweep interval, period, and perigee change) need adjustment. 

* Of course the analysis may be rendered invalid if the station is in the orbit 
plane at the time of any of the 3 data points. The first point was reported only as a 
matter of convenience. In any event the comparison of predicted and measured 
trajectories described later serves as a final check on the validity of the elements. 



ORBITAL ELEMENT GENERATION 647 

Fig. 14 shows the fit of an orbit generated from July 30 data only 
(azimuth, elevation, range). Here the pointing fit is within 0.02° for that 
day, but degrades to 0.4° a day later and is grossly inadequate when 
matched with data from the June 2 trajectories. Again better rate in- 
formation is needed. 

In Fig. 15, we employ the MOERATE program output to determine 
the orbital rates needed to connect the June 30 and July 30 elements 
(see Sections IX and 10.1.2). With but a few comments this output should 
be self-explanatory. First of all, the theoretical values at the bottom of 
the page are computed from the set of MOES listed at the top left of the 
figure. The chief purpose for their computation is to gain an estimate of 
the rates for the actual rate calculations. One cannot compute an 
average anomalistic period from two given epochs if an estimate of the 
period is not available. Without such an estimate, the satellite may have 
made any number of passes through perigee between the two epochs and 
an infinite array of average periods would be valid solutions. Even with 
such an estimate, it is sometimes possible to miscalculate the number of 
times a satellite has passes through perigee if the time between epochs 
becomes sufficiently large. That is why the number of perigee passes 
computed by the program is printed out. Fig. 15 (see checked line) has 
been chosen to show a misestimate of one in the calculation of the peri- 
gee passages (= 196) and the attendant error in the average period. By 
increasing the perigee passes to 197 and performing the indicated divi- 
sion, a new average period of 225.30083 is calculated which is closer to 
that reported by the first set of MOES and is, of course, a valid solution. 
Strictly speaking, if this new period value is used, the apsidal advance 
per T a should be corrected by multiplying it by 225.30083/226.45032 to 
obtain 0.19129°/T a . The PSI is unaffected by these operations. 

Printed at the very bottom of Fig. 15 is the number of iterations re- 
quired for the computed PSI and apsidal advance rate to come within 
2.1 minutes and 0.2 degrees of the theoretical rates. The program stops 
after 21 iterations whether or not the convergence is successful. 

Since there is no reason to suspect that the anomalistic period associ- 
ated with the June 30 orbit generation possesses high validity (see Fig. 
13), the portion of the MOERATE output dealing with changes in period 
is ignored and the indicated average anomalistic period is used. From Fig. 
15, the computed rates become 

anomalistic period = 225.30083 minutes 

apsidal advance rate = 0.19129° per anomalistic period 

prime sweep interval = one day — 8.12511 minutes. 

These rates are then submitted along with the June 30 data to MOE- 
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ORBITAL ELEMENT GENERATION G53 

GEN G3. The output is shown in Fig. 16. A considerable improvement 
in the fit over the June 2 to August 1 period is seen, with the bulk of the 
errors within the low hundredths of degree range. A maximum range 
error of only 7 miles over this interval is also noted. 

Even these small errors may indicate the need for additional orbital 
rate improvement. The mean rates determined by combining June 2, 10, 
30 and July 30 data are 

anomalistic period = 225.30049 minutes 

apsidal advance rate = 0.19105° per anomalistic period 

prime sweep interval = one day — 8.12480 minutes. 

Using these rates produces an even better trajectory fit, as shown in Fig. 
17. The mere fact that the maximum measured error is about 0.05° and 
within five miles in range over a 3-month period indicates the near con- 
stancy of the secular rates and attests to the fact that little energy is 
being lost by the satellite as it orbits. 

Fig. 18 illustrates the end result of the same type of procedure as de- 
scribed above, except that only one range point per determination is 
used in the generation. Here the pointing errors peak at 0.6° over the 
survey with maximum range errors about 25 miles.* 

Fig. 19 repeats the above with absolutely no range data fed to the 
orbital generation process. Over the survey, maximum pointing and 
range errors are 0.8° and 76 miles respectively. 

The advantage of 3 range points per determination over the several 
months' survey period is clearly evident by the results. 

It must be stated that this illustrative example is somewhat simplified 
in that a certain amount of initial, free-fit orbit generation is omitted. 
In actual operation, the free-fit mode of MOEGEN G3 is often run for 
6 to 12 data points in various combinations, and only the elements pro- 
ducing the best trajectory fit are preserved for later rate processing. 

At this point, we intend to show how the techniques of Section VIII 
will provide a better round-the-world orbit fit by modifying the eccen- 
tricity and perigee radius of the orbital elements. Since no data pertain- 
ing to the TELSTAR satellites were readily obtainable from stations in 
southern latitudes, it was decided to proceed as described below. 

A set of orbital elements f (which we shall call MO) was generated 



* The -7449.740 and -7211.212 range "errors" indicated by Fig. 18 are not 
true errors but merely reflect the fact that no range data were supplied at these 
times for generation purposes. The true errors are ±3 miles, as shown several lines 
below. 

t These orbital elements had previously predicted Andover passes to a pointing 
accuracy of 0.07°. 
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Table II — Alterations 


to MO Elements 


Data Point of 
4th Trajectory 


Azimuth 
Alteration 


Elevation 
Alteration 


Range 
Alteration 


1 
2 
3 


0.21° 

0.06° 

-0.10° 


-0.07° 
-0.02° 
+0.02° 


4.85 miles 

1.03 miles 

-0.87 mile 



from Andover data and used to simulate trajectories over Andover for a 
sequence of four passes. The computed pointing angles for these passes 
then were rounded off to the nearest hundredth of a degree so as to re- 
semble typically reported data from Andover. Three points from each 
of the first 3 passes were then selected for submission to the MOEGEN 
G3 program. Three points from the fourth pass were purposely altered 
before submission to simulate the effects of tracking errors, mechanical 
jitter of the ground antenna, errors in assessment of refraction which 
affect measured position, and digital noise in the equipment at the 
ground station. These alterations were slight and are given in Table II. 
Using the data from the 4 passes, MOEGEN G3 produced 4 sets of or- 
bital elements. 

Fig. 20 displays the generated elements Ml, M2, M3 and M4. For 
each of these elements, the first 3 lines of data under "pointing angle 
comparisons" (see Fig. 20) show the differences between the input data 
points of the submitted trajectories and those predicted by the derived 
elements. For each set of elements, the pointing error is within 0.07° and 
the range difference under 0.9 miles. Looking at each of the first 3 com- 
parison lines, all 4 sets of elements produce excellent fits to reported 
Andover trajectories, and from this examination there would be no rea- 
son to expect that any particular set of elements would be decidedly 
better than any other. However, comparisons of the M4 orbit with 
6/30 trajectory data indicate that something is wrong with these ele- 
ments because of the increase in pointing angle differences on this date. 
This, however, could be produced by variety of causes: e.g., errors in 
epoch, rates, ellipse shape, size and orientation of the orbit. The point to 
make is that the M4 orbit fits well the trajectory over Andover from 
which it was derived; the purposely introduced errors in that trajectory 
have not done violence to the orbit theoiy, so that the motion of the 
satellite is in essential accord with those motions permitted by the earth's 
gravitational potential. 

Using the same MO elements as previously, trajectories over Johannes- 
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burg, South Africa were produced. These were compared with the pre- 
dictions from the generated M4 elements in a computer comparison pro- 
gram. The results are given in Fig. 21, where all listed data are for 
satellite positions not visible to Andover (subsatellite latitudes 4° to 
—42°). The average pointing error is 2.7° and the root-mean-square 
range error is 65.5 miles. 

The four generated elements were then submitted to the MOEFIT 
program, which applies the principles of Section VIII to generate an 
ellipse of best fit using only those portions of the orbits, described by the 
input elements, which are visible from Andover. Fig. 22 shows the ec- 
centricity of 0.400593 and the perigee radius of 4565.168 statute miles 
determined from this ellipse. It also indicates that the root-mean-square 
departure of the fitted ellipse from those points generated by the four 
elements is in the order of 10 -7 . Fig. 23 graphs the fitted ellipse as a curved 
line, while the points visible from Andover and predicted by the in- 
dividual elements are noted by zeros on the plot. 

The M4 set of elements was then altered by substituting the eccen- 
tricity and radius of perigee values from the MOEFIT program and 
was again compared with Johannesburg trajectories given by the original 
set of elements, MO. Typical results are given in Fig. 24, where it is 
seen that the average pointing angle error is now 0.07°, about 2 orders 



RIUND-THE-WBRLD-BRB1T-TEST - MBES GENERATED FM AND0VER L C THflHAS 

BRBMAL ELEMENTS N-4 FRIN HBEGEN G3 USED. 

STATIBN AT JBBURG, S.AF LAT- -26. 03246 DEG WEST LBNG - 331.75878 DEG HI AB0VE HSL= 5370.00000 FT 



YEAR 


MB OV 


HR 


MN 


SEC 


SAT PBS FRBf 


1 MBE 1 SAT 


PBS FRBM DATA 


D I F 


FERE 


N C E S 














AZ-CEG 


EL-DEG 


RANGE-HI AZ-OEG 


EL-DEG 


RANGE-MI 


AZ-DIFF 


EL-DIFF 


ARC-DIF RANGE-OIFF 


1964 


6 30 


t 


50 





288.75 


11.10 


7605.79 288.53 


11.98 


7559.01 


-0.221 


0.878 


0.9043 


-46.777 


1964 


6 30 


2 








284.15 


18.01 


6463.13 283.78 


19.20 


6411.44 


-0.367 


1.191 


1.2407 


-51.693 


1964 


6 30 


2 


10 





277.26 


26. 86 


5135.54 276.55 


28.50 


5085.55 


-0.713 


1.635 


1.7529 


-49.993 


1964 


6 30 


2 


20 





263.52 


39.54 


3647.60 261.69 


41.83 


3616.87 


-1.829 


2.287 


2.6741 


-30.732 


1964 


6 30 


2 


30 





213.29 


54.33 


2253.39 206.47 


55.47 


2295.77 


-6.816 


1.137 


4.0786 


42.382 


1964 


6 30 


2 


40 





136.40 


14.84 


2296.01 134.91 


14.78 


2430.20 


-1.487 


-0.057 


1.4384 


134.186 


1964 


7 1 


4 


40 





251.49 


2.50 


5127.70 250.75 


4.33 


5068.46 


-0.737 


1.834 


1.9758 


-59.242 


1964 


7 1 


4 


50 





236.65 


14.72 


3173.39 235.10 


17.25 


3136.97 


-1.549 


2.525 


2.9315 


-36.423 


1964 


7 1 


5 








174.13 


27.48 


1518.65 169. 84 


29.93 


1586.01 


-4.293 


2.448 


4.4906 


67.362 


1964 


7 1 


a 


56 





277.23 


31.34 


832.02 277.78 


37.04 


906.03 


0.549 


5.697 


5.7154 


74.006 


HCDIF1EC BR 


BITAL 


ELEMENTS FBR 


ABBVE 


TABLE ARE- 








AVG 


=2.7202 


6S.S4I=RMS 


CLEMENTS NB 


1 






















REFERENCE TIME 1964 


7 1 5 


HR 11. 


13300 MIN UT 















INCLINAJIBN 42.72288 DEG 
ASCENOING NBOE WEST LBNG 255.82689 DEG 
PRIME SWEEP INTERVAL ' 1 DAY -8.124800 MIN 
PFRIGEE AND SATELLITE ARGUMENT 324.12966 DEC 

RATE BF CHANGE 0.191050 DEG PER PERIBD 
ANCMALISTIC PERIBD 225.30049 MIN 

RATE BF CHANGEO. MIN PER PERIBD 

ECCENTRICITY 0.4115B10 
RADIUS IF PERIGEE 4450.926 STATUTE MILES 



Fig. 21 — Comparison of M4 elements to Johannesburg trajectories derived 
from MO. 
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R0UND-THE-W0RLO-0RBIT-TEST - M0ES GENERATEO FM ANDOVER L C TH0MAS 
SUCCE5SI0N OF PASSES 0N 6/30 ANO 7/1/1964 USEO F0R RP.E DETERM1NATI 0N 

ABS VALUE 0F MAX ERROR- 5 .8 8 22". 79E-0 7 
LEAST SQUARES ORBITAL ELEMENT VALUES- MEAN SQUARE ERROR* 2 . 1384752E-14 

RMS ERROR- 1.4623526E-07 
ECCENTRICITY « 0.400593 RADIUS 0F PERIGEE = 4565.1683 STATUTE MILES 
SEMI-MAJ0R AXIS = 7616.1472 STATUTE MILES 

AVERAGE PERIGEE POSITION = 324.14983 DEC AT EP0CH 0F FIRST M0ES 

INPUT M0OIFIEO 0RBITAL LIST FOLL0WS- 

REFERENCE TIME 1964 6 30 17 H0UR 55.21*00 MIN U.T. 

INCLINATI0N 42.74526 OEG 

ASCENOING N0DE L0NG 85.920820EG WEST 

PRIME SWEEP INTERVAL 0NE OAY -8.12*80 MIN 

PERIGEE ANO SATELLITE P0SITI0N 323.58760 OEG 

RATE 0F CHANGE 0.19105 OEG PER PERIOD 
ANBMALISTIC PERIOD 225.30049 MIN 

RATE 0F CHANGE 0. MIN PER PERI0D 
ECCENTRICITY 0.40099 
RADIUS 0F PERIGEE 4564.55493 STATUTE MILES 



START CARD IMAGE I S- 
1964 6 30 18 7 



70 60 60 



REFERENCE TIME 1964 6 30 21 H0UR 40.50300 MIN U.T. 

1NCLINATI0N 42.74677 DEG 

ASCENDING N0DE L0NG 142.55847DEG WEST 

PRIME SWEEP INTERVAL 0NE OAY -8.12480 MIN 

PERIGEE AND SATELLITE P0SITI0N 323.75162 DEG 

RATE 0F CHANGE 0.19105 OEG PER PERI0D 
ANOMALISTIC PERI0O 225.30049 MIN 

RATE 0F CHANGE 0. MIN PER PERIBD 
ECCENTRICITY 0.40110 
RADIUS 0F PERIGEE 4563.15997 STATUTE MILES 



START CARD IMAGE I S- 
1964 6 30 22 1 



124 120 120 



H0UR 25.83200 MIN U.T. 



REFERENCE TIME 1964 7 1 

INCLINATION 42.74387 OEG 

ASCENOING NODE LONG 199.19972DEG WEST 

PRIME SWEEP INTERVAL ONE DAY -8.12480 MIN 

PERIGEE AND SATELLITE POSITION 323.98413 DEG 

RATE OF CHANGE 0.19105 DEG PER PERIOD 
ANOMALISTIC PERIOD 225.30049 MIN 

RATE OF CHANGE 0. MIN PER PERIBD 
ECCENTRICITY 0.40093 
RADIUS 0F PERIGEE 4566.15997 STATUTE MILES 



Iffl 



START CARD IMAGE I S- 
1964 7 12 2 



138 120 120 



REFERENCE TIME 1964 7 1 5 HOUR 11.13300 MIN U.T. 

INCLINATION 42.72288 DEG 

ASCENDING NODE LONG 255.826B9DEG WEST 

PRIME SWEEP INTERVAL ONE OAY -8.12480 MIN 

PERIGEE ANO SATELLITE POSITION 324.12966 OEG 

RATE OF CHANGE 0.19105 DEG PER PERIOD 
ANOMALISTIC PERIOD 225.30049 MIN 

RATE 0F CHANGE 0. MIN PER PERIOD 
ECCENTRICITY 0.41158 
RADIUS 0F PERIGEE 4450.92596 STATUTE MILES 



CO 



START CARD IMAGE I S- 
1964 7 1 6 14 



99 60 60 



Fig. 22 — Eccentricity and perigee radius from ellipse of least-squares fit. 



of magnitude better than the value before ellipse fitting. The root- 
mean-square range error has decreased an order of magnitude and is now 
5.2 statute miles. 

It is to be noted that only four passes at Andover were used in this il- 
lustrative procedure, since these four encompassed the entire portion of 
the orbit visible from Andover. However, as apsidal advance takes place, 
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new portions of the orbit become visible and, provided little orbit decay 
occurs, these may also be used to obtain even better ellipse fits and fur- 
ther refine values for the eccentricity and perigee radius. 

It is also noted that in actual practice each of the trajectories used 
would have suffered slight measurement errors akin to those artificially 
introduced into the fourth trajectory that produced elements M4. None- 
theless, the derived elements would predict well that portion of the orbit 
ellipse from which they came (as was demonstrated by M4) and would 
therefore be improved by MOEFIT for round-the-world use by about 
the same amounts as shown in the illustrated example. 
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Fig. 23 — Ellipse of least-squares fit to elements. 
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XII. CONCLUSIONS 

The methods described have been used extensively for those orbit 
determinations of the TELSTAR satellites dealing with the satellite 
communications experiment. The following merits of this orbital tech- 
nique are noted: 

(i) The elements chosen explicitly express the secular rates, thereby 
simplifying and making more economical the generation of predicted 
trajectories and drive tapes for the Andover antenna. 

(ii) The generation procedure is geared to sparse and discontinuous 
tracking data. 

(in) Using a modified Gaussian technique for orbit generation pro- 
vides a straightforward method for range determination to which orbital 
perturbations are easily added on an iterative basis. 

(iv) The computation of secular rates only saves time and money and 
proves adequate for the tracking requirements of the TELSTAR satellite 
experiments (0.05° pointing accuracy). 

(v) For cases where the secular rates are not constant over a given 
span of time, the methods described will still function provided tracking 
data are available frequently enough to assure essentially constant 
rates from each data acquisition time to the next. 

(vi) The method provides for the operational inclusion of perturbations 



RIUND-THE-WBRLD-IRB1T-TEST - MBES GENERATED FM ANDBVER I C THBMAS 

M-4 ELEMCNTS WITH ECCENTRICITY AND RADIUS 0F PERIGEE FR6M NBEFIT. 

STATIBN AT J B BURG, S.AF LAT> -26.C3246 DEG WEST LBNG • 331. 75878 DEG HT ABBVE MSL- 5370.00000 FT 



YEAR MB 0Y HR MN SEC 



1964 
1964 
1964 
1964 
1964 
1964 
1964 
1964 
1964 
1964 



6 30 1 50 

6 30 2 

6 30 2 10 

6 30 2 20 

6 30 2 30 

6 30 2 40 



SAT PB 
A7.-CEG E 
2BB.56 
283. 81 
276.59 
261.72 
206.38 
134.85 
250.76 
235.13 
169.84 
277.52 



S FRBH HBE 
L--DEG RANGE 
11.97 7550 



19.20 
28.51 
41.85 
55.52 
14.72 
4.16 
17.22 
29.94 
36.95 



6403 
5078 
3611 
2292 
2429 
5062 
3133 
1584 
906 



1 SAT 
-MI AZ-0EG 
,87 2B8.53 
,74 283.78 
,57 276.55 
,05 261.69 
,20 206.47 
,98 134.91 
,97 250.75 
,32 235.10 
.84 169.84 
.26 277.78 



PCS FRBM DATA 

EL-DEG RANGE-MI 

11.96 7559.01 



I F F E R E 

A2-DIFF EL-DIFF 

-0.027 0.009 



N C E S 

ARC-DIF RANGE-D1FF 
0.0280 8.142 



19.20 
28.50 
41.83 
55.47 
14.78 
4.33 
17.25 
29.93 
37.04 



NIDIFIED BRBITAL ELEMENTS FBR ABBVE TABLE ARE- 
ELEMENTS NB 1 
REFERENCE TIME 1964 7 1 5 HR 11.13300 MIN UT 
INCLINATIBN 42.72288 DEG 
ASCENDING NBDE WEST LBNG 255.B26B9 OEG 
PRIME SWEEP INTERVAL • 1 DAY -8.124800 MIN 
PERIGEE AND SATELLITE ARGUMENT 324.12966 DEG 

RATE BF CHANGE 0.1910SO OEG PER PERISO 
ANBMALISTIC PERIBO 225.30049 MIN 

RATE BF CHANGEO. MIN PER PERIBD 

ECCENTRICITY 0.4005930 
RADIUS BF PERIGEE 4565.168 STATUTE MILES 



6411.44 
5085.55 
3616.87 
2295.77 
2430.20 
5068.46 
3136.97 
1586.01 
906.03 



-0.027 

-0.040 

-0.034 

0.091 

0.055 

-0.029 

-0.026 

0.000 

0.264 



0.003 

-0.007 

-0.017 

-0.047 

0.064 

0.168 

0.026 

-0.012 

0.093 



0.0242 
0.0370 
0.0321 
0.0699 
0.0833 
0.1700 
0.0350 
0.0171 
0.2305 



7.704 
6. 983 
5.821 
3.566 
0.222 
5.492 
3.653 
1.166 
-2.234 



A VG =0. 0727 S. 209 =RMS 



Fig. 24 — Comparison of M4 elements (with fitted ellipse corrections) to 
Johannesburg trajectories. 
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not computed here. For example, cyclic perturbations may be expressed, 
if indicated in the data, by generation of a multiplicity of elements over 
a time span and noting the derivatives of the "secular" rates. 
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APPENDIX A 

Relationship of the Triangular Area between Two Radius Vectors and the 
Orbit Normal 

Consider triangle 023 of Fig. 4 and project it upon the yz plane. The 
vertices of the projected triangle are, therefore, (0,0,0), (0,7/ 3 ,z 3 ), 
(0,2/2,22) and hence its area becomes 

A-v = M2/2Z3 - Z22/3). (119)* 

Now the normal to the yz plane is, of course, the x axis. The normal to 
the orbit plane, shown in Fig. 4, makes an angle n with the x axis. Since 
the area of triangle 023 is the projected area referred to above divided 
by the cosine of the angle between the orbit plane and the yz plane, and 
since this angle is the same as the angle n between the normals to the 
two planes, it follows directly that the area of 023 is 

A m = M2/2Z3 - z 2 y 3 )/cosn (120) 

or 

2/2Z3 — z-iy-A = A Ti cos n (9) 

where 

^23 = 2^023- (121) 

APPENDIX B 

Relations Involving the Earth's Oblateness 

Because the earth is not a perfect sphere, two corrections are normally 
applied to satellite positional data as measured from any specified ground 



* See page 266, Ref. 23. 
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station. The same corrections in the inverse sense are applied when tra- 
jectory predictions are made from orbital elements. The first of these 
recognizes that the geocentric distance of the station is a function of 
geodetic latitude; the second correction relates to the change in azimuth 
and elevation of the satellite when referred to a spherical as opposed to 
an oblate earth. 

The geocentric distance 5 of the station is simply 

R = H + R E (0.998320047 + 0.001683494 cos 2<p 

(125) 
- 0.000003549 cos 4*> + 0.000000008 cos 6^) 

where 

H = the height of the station above mean sea level in statute 

miles 
Re = the earth's equatorial radius = 3963.347 statute miles 
<p = the geodetic latitude of the station. 

This indicated value of R is used in solving for the satellite's slant range 
as in Section IV of the text. It is also used to determine the X, , Y ',• , Z, 
station coordinates in Section 5.3. 

Normally, azimuth-elevation measurements of a satellite made from 
a ground station are based upon measurements referenced to a local 
horizontal plane. This plane is often taken as a tangent to the oblate 
earth at the station's location. It becomes useful to reference such azi- 
muth-elevation measurements to a plane tangent at the station to a 
sphere having a radius R. This is done both in trajectory prediction (Sec- 
tion IV) and in preprocessing the azimuth-elevation data prior to orbit 
determinations. 

To relate the azimuth-elevation points to the two tangent planes, we 
begin by establishing two right-hand coordinate systems, centered on 
the station. In the first of these, the x-y axes lie hi the plane tangent to 
the oblate earth with the y axis pointing northward and the x axis east- 
ward. In the second system, the x-y axes lie in the plane tangent to a 
sphere of radius R. The transformation from one system to the other is 
a simple one, since rotation of the first set of coordinates about the com- 
mon x axis moves y into y as shown in Fig. 25(a). 

By inspection of Fig. 25(b), we write the measured azimuth (A) and 
elevation (E) in terms of the xyz coordinates (referenced to the oblate 
earth) as 

x = p cos E sin A (126) 

y = p cos E cos A (127) 

z = p sin E (128) 
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y »' 



NOTE: 
L) IS TANGENT TO 
OBLATE EARTH 

L)' IS TANGENT TO 
SPHERE OF 
RADIUS R 



P- STAT I ON 



(a) 



(b) 




P COS E SIN A 



where 



Fig. 25 — Geometry for correcting earth oblateness effects. 



p = the slant range to the satellite. 



Next we rotate the xyz coordinates about x through an angle (a)* equal 
to the difference between the geodetic and geocentric latitude of the 
station. This brings y into y and z into z where xy z is the coordinate 

* a. = geocentric latitude — geodetic latitude « 0.1962 sin (geodetic latitude). 
A closer approximation is found in Ref. 5. 
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system for a sphere of radius R. The relationships between y,y and z,z 
are 

y = y cos a + z sin a (129) 

z = —y sin a -\- z cos a. (130) 

Combining (126) through (130) one obtains 

x = p cos E sin A (131) 

y = p cos i£ cos A cos a + p sin i£ sin a (132) 

z = — p cos E cos A sin a -\- p sin 2? cos a (133) 

from which the azimuth (A') and elevation (E') referred to the sphere 
become 

E' = sin -1 zip (134) 

A' = tan -1 x/y'. (135) 

The procedure is perfectly symmetrical, and if one wants to convert 
from spherically oriented azimuth-elevation to the oblate earth azimuth- 
elevation, (131) through (135) are used with a negative a and inter- 
changed primes. 

appendix c 

The Triangle Ratios nil and m 3 as Time Functions 

Consider the orbit ellipse lying in the xy plane as shown in Fig. 26. 
The observed positions of the satellite are expressed in x,y coordinates and 
time /. If we now define a new time variable t = kt, the familiar differ- 
ential equations of motion then become 



and 



!&=-*' 3- (139) 

dr- r 



Considering points 1 and 2, we may expand Xi and yi in a Taylor series 
which generally converges for all reasonable values of t: 

a ; 1 = ^ 2 + $ T3 + i-^ T 3 2 +... (140) 



dr ' ' 2\di 
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y 



Gc 3 >y 3 >t. 3 ) 



FOCUS 



(«2,y2,t 2 ) 



(xi,y,,t,) 




Fig. 26 — Orbit triangles in x,y,t space. 

dij 1 d 2 y 2 , 

1Jl = V2 + d-T Ts + 2\d? T3 + '"' 



(141) 



where the r derivatives are to be evaluated at the time t 2 of the second 
data point. We define t 3 in the text. Similar expressions may of course be 
written for points 2 and 3. It is possible to substitute for all the deriva- 
tives in (140) and (141) higher than the first by successive differentia- 
tion of (138) and (139). This permits (140) and (141) to be expressed 



as 



xi = Aix 2 + Bi -j- 



yx = A x y 2 + B t 



dy 
dr 



where 



Ai = 1 - ^ 



2 \r 2 3 r 2 4 dr) 



+ 



7-3 



T3 



dr 



B\ — tz — n — \ t i — ; t "i 
6r 2 3 4r 2 4 dr 



(142) 



(143) 



(144) 



(145) 



Point x 3 , 2/3 may be expressed in a manner identical to (142)- (145) with 
subscript 1 being replaced by 3 throughout and t 3 becoming r 2 . 
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The areas of the triangles formed by n , r 2 , and r 3 and their respective 
chords are 

A 12 = \(xjjf% - 2/1*2) (146) 

and 

An = Hwh - ViXz). (147) 

If now (142) and (143) are substituted into (146) then, 

*.-»(*£-»!) 048) 

But the right-hand member of (148) is B t times the areal velocity* of 
the satellite, which by Kepler's Third Law is simply 

A u - §Bifc Vp (149) 

where 

p = a(l — e 2 ) 

fc = the Gaussian constant for 
earth satellites 

= 0.07436574 min" 1 . 

Similar expressions exist for An which parallel (148) and (149) with sub- 
script 3 substituting for 1. The area enclosed by n , r 3 and the connecting 
chord differs a bit from either A 12 or A23 and is 

An - £(*i2/3 - 2/1*3) 

_ (150) 

= \{A,B, - BxAJky/p- 

Now by explicitly expressing the A { , J5, factors of (144) and (145) in 
the sector area equations (149) (150), we have 

^--»Vr.[i-^ + ^fe^J+-]- 053) 



See Ref. 14; Ref. 24, p. 5; and Ref. 9, p. 26. 
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The mi and m 3 equations follow directly from the above as 

Mi = -j- 

_ ti I" , t 3 (t 2 + n) t 3 (t3 2 + tit 3 — ti 2 ) dr , ~| 



m 3 = 



A u 

A 13 



(155) 

_ T3 T Ti(t 2 + T 3 ) _ Ti(ti + T1T3 — T3") rfr , I 

" r 2 |_ + 6r 2 3 4r 2 < dr ~ l ~ "J* 
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